//Copyright 2010 Neurosciences Research Foundation, Incorporated
/*
 *  microcircuitry.h
 *  
 *
 *  Created by Jeff McKinstry on 6/26/09.
 *
 *  NOTE: microcircuitry depends completely on neuronspecs.dat file.  If the neurons aren't there, you can't connect them.
 *  Also, layer info is still hardcoded in neuronspec.h for now.  
 */
#ifndef MICROCIRCUITRY_H
#define MICROCIRCUITRY_H

#include <vector>
#include <string>
#include <iostream>
#include <fstream>
#include "neuronspecs.h"
#include "areas.h"
using namespace std;



//**************************************************************
/*
 Reads microcircuitry.dat file to set up the local connectivity.
 Pre: read_groups() has been called.
 Post: SYN table and synapses table have been filled.
 
 // distribution of synapses for different layers for different neuronal types
 int	synapses[Ntypes][Layers]={}; // Read in from microcircuitry.dat
 
 vector<float>	SYN[MAX_AREAS][Ntypes][Layers];  
 // To save space, all vectors of synaptic contact percentages are zero, unless microcircuitry.dat has connections specified in that area.
 
 string layerNames[] =  {"L6", "L5", "L4", "L23", "L1", "Th", "RTN"} is used to translate layer names to layer numbers.
 
 
 Anatomy matrix. Cortico-cortical connections, brainstem, and sensory connections are in the last column.  The last two are for future use.
 The file that fills this structure takes on the following format, similar to the table in Izhikevich and Edelman, 2008 PNAS paper, supporting information.
 c-c   brstm  sensory
 nb1	L1		8890	10.1,  6.3,  0.6,  1.1,  0.0,  0.0,  0.1,  0.0,  0.0,  0.1,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  4.1,  0.0,  0.0,  0.0, 77.6,  0.0,  0.0, 
 p23	L23		5800	 0.0, 59.9,  9.1,  4.4,  0.6,  6.9,  7.7,  0.0,  0.8,  7.4,  0.0,  0.0,  0.0,  2.3,  0.0,  0.0,  0.8,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0,
 ...
 
 
 If min_radius_mm is zero, then the connection will be a surround, donut-shaped connection with a Gaussian centered midway between min and max radius.
 
 */


// ST-STDP types:
enum STSTDP_type {NONE, REGULAR, LINEAR_REDUCTION, POLY_REDUCTION};


class connection  {
public:

	connection(){
		pct_pre_synapses = 0.0;
		sigma_mm = 0.0;
		min_radius_mm = 0.0;
		max_radius_mm = 0.0;
		init_noise_pct = 0.0;
		init_max_sum = 0.0;
		max = 0.0;
		lrate0 = 0.0;
		lrateN = 0.0;
		lstarttimems = 0;
		lendtimems = 0;
		DAstarttimems = 0; //RGM: Added on 02/22/2012
		DAendtimems = 0;
		ampa_gain = 0.0;
		nmda_gain = 0.0;
		nmdaVi_gain = 0.0;
		gabab_gain = 0.0;
		ststdp_type = NONE;
		ststdp_modulation = 0.0;
		ststdp_max = 0.0;
	}

	
	float pct_pre_synapses;
	float sigma_mm;
	float min_radius_mm;
	float max_radius_mm;
	float init_noise_pct;
	float init_max_sum;
	float max;
	float lrate0;
	float lrateN;
	int lstarttimems;
	int lendtimems;
	int DAstarttimems; //RGM: Added on 02/22/2012
	int DAendtimems;
	float ampa_gain;
	float nmda_gain;
	float nmdaVi_gain;
	float gabab_gain;
	int ststdp_type;
	float ststdp_modulation;
	float ststdp_max;
	
};


//*********************************************
template<typename T>
int find(T key, T array[], int max){
	
	int index = -1;
	for(int i = 0; i < max; i++){
		if( array[i] == key ){
			index = i;
			break;
		}
	}
	return index;
	
	
}





template <int MAX_AREAS, int MAX_TYPES, int Layers>
class microcircuitry{
	
public:
	
	
	microcircuitry(){
		
		neuronspecs ns;
		layerNames = ns.get_layer_names();
		ns.read_neuronspecs();
		Names = ns.get_neuron_names();
		Ntypes = ns.get_num_types();

		
		for(int area = 0; area < MAX_AREAS; area ++ ){
			for(int i = 0; i < MAX_TYPES; i++){
				for(int j = 0; j < Layers; j++){
					synapses[area][i][j] =0 ;
				}
				cell_type_percents[area][i]=0.0f;

			}
		}
		
		
	}
	
	// pre: read_microcircuitry has been called.
	// post: synapses and SYN structures copied.
	void	get_microcircuitry(int	out_synapses[MAX_AREAS][MAX_TYPES][Layers], vector<vector<connection> >	out_SYN[MAX_AREAS][MAX_TYPES][Layers]){
		
		for(int area = 0; area < area_data.get_N_areas(); area ++ ){
			for(int i = 0; i < Ntypes; i++){
				for(int j = 0; j < Layers; j++){
					out_synapses[area][i][j]= synapses[area][i][j];
				}
			}
		}
		
		for(int area = 0; area < area_data.get_N_areas(); area ++ ){
			for(int i = 0; i < Ntypes; i++){
				for(int j = 0; j < Layers; j++){
					if( SYN[area][i][j].size() > 0 ){
						out_SYN[area][i][j] = SYN[area][i][j];
					}
				}
			}
		}
		
	}
	
	void get_neuron_percentages(float out_cell_type_percents[], string area_name){
		
		int area = get_area_number(area_name);
		
		if(area < 0 ){
			cout << "Class microcircuitry: function get_neuron_percentages: area name " << area_name << " was not found in microcircuitry.dat." << endl;
			exit(1);
		}
		
		for(int i = 0; i < Ntypes; i++){
			out_cell_type_percents[i] = cell_type_percents[area][i];
		}
		
	}
	
	// returns -1 if area_name not found.
	int get_area_number(string area_name){
	
		return area_data.get_area_number(area_name);
		
	}
	
	//************************************************
	void read_microcircuitry(int myid=0){
		
		
		area_data.load_areas("areas.dat");
		
		ifstream fin("microcircuitry.dat");
		
		string area_name, first_string, cell_type_name, layer_name;
		int num_synapses, layer;
		float cell_type_percent;
		
		// Read the area name. It must exist in areas.dat file also.
				
		int i = 0;
		if (!fin.is_open()) cerr << "OOOPS, didn't open microcircuitry.dat!" << endl;
		while ( fin >> area_name ){
			if ( area_name[0] == '%' ) // this is a comment line, skip it
			{
				cout << "discarding comment" << endl;
				string foo;
				getline( fin, foo); // discards the rest of the line
				// get beginning of next line
				continue;
			}

			if( ! area_data.exists(area_name) ){
				cout << "Error reading microcircuitry.dat areas. " << area_name << " is an unknown area name @ char " << fin.tellg() <<".  It was not found in 'areas.dat'\n";
				exit(0);
			}
			
			cerr << area_name << endl;
			int area = area_data.get_area_number(area_name);
			

			fin >> first_string;
			
			int type;
			while( first_string[0] != '*' && !fin.eof()){
			  

				if ( first_string[0] == '%' ) // this is a comment line, skip it
				{
				  cerr << "discarding comment" << endl;
				  string foo;
				  getline( fin, foo); // discards the rest of the line
				  fin >> first_string; // get beginning of next line
				  continue;
				}

				if( first_string[0] != '"' ){
					// First line for this post-type.  Read in the cell parameters
					cell_type_name = first_string;
					
					fin >> cell_type_percent;
					
				
					//cout << "cell_type_name " << cell_type_name[0] << " " << layer_name << " " << cell_type_percent << " " << num_synapses << "; ";
					// lookup the type number
					type = find( cell_type_name, Names);
					if( type < 0 ){
						cout << "Error reading microcircuitry.dat column 1. " << cell_type_name << " is an unknown cell type name."<< endl;
						cout << "Line: " << cell_type_name << " " << layer_name << endl;
						exit(0);
					}
					//cout << "type " << type << " ";

					// store cell_type_percent.
					if( cell_type_percent > 0.00000001 ) cell_type_percents[area][type] = cell_type_percent;  // The first entry for a given cell type must contain the percentage of
				}
				else{
					// Subsequent lines for this post-type have " marks in the first 2 spots.  Get rid of the second spot.
					string foo;
					fin>> foo;
				}

				// Now read in the next 2 fields, the layer and the #synapses in that layer.  Either they are both quotes, or they are both defined.
				fin >> layer_name;
				if( layer_name[0]!='"'){
					fin >> num_synapses;

					// lookup the layer number
					layer = find( layer_name, layerNames);
					if( layer < 0 ){
						cout << "Error reading microcircuitry.dat column 2. " << layer_name << " is an unknown layer name."<< endl;
						cout << "Line: " << cell_type_name << " " << layer_name << endl;
						exit(0);
					}
					//cout << "layer " << layer << " " ;
					
					// this type of cell for this area;  all subsequent lines for this cell type in a given area must be set to zero.
					if (myid == 0)cerr << num_synapses << " " ;
					synapses[area][type][layer] = num_synapses;
					
					//if(myid==0)cerr << "cell_type_percent: " << cell_type_percent << endl;
					
					// resize
					SYN[area][type][layer] = vector<vector<connection> >(MAX_AREAS, vector<connection>(Ntypes));
					
				}
				else{
					// Subsequent lines for this post-type have " marks in the first 2 spots.  Get rid of the second spot.
					string foo;
					fin>> foo;
				}
								

				// Read in the next 2 fields: pre_area_name and pre_type_name and translate to pre_area and pre_type numbers for array indexing.

				string pre_area_name, pre_cell_type_name;
				fin >> pre_area_name >> pre_cell_type_name ;
				
				int pre_area;
				
				if( pre_area_name[0] == '-'){
					// Default is indicated by a dash.  Change to the current area.
					pre_area = area;
				}
				else{
					if( ! area_data.exists(pre_area_name) ){
						cout << "Error reading microcircuitry.dat pre areas. " << pre_area_name << " is an unknown area name @ char " << fin.tellg() <<".  It was not found in 'areas.dat'\n";
						exit(0);
					}
					pre_area = area_data.get_area_number(pre_area_name);
				}
				
				// lookup the type number
				int pre_type = find( pre_cell_type_name, Names);
				if( pre_type < 0 ){
					cout << "Error reading microcircuitry.dat column 6. " << pre_cell_type_name << " is an unknown cell type name."<< endl;
					cout << "Line: " << cell_type_name << " " << layer_name << endl;
					exit(0);
				}
				//cout << "type " << type << " ";
				
				

				float pre_percent, min_radius_mm, max_radius_mm , sigma_mm, init_noise_pct, init_max_sum , max, lrate0,lrateN;//RGM: added tau_dopamine on 02/22/2012
				int lstarttimems,lendtimems, DAstarttimems, DAendtimems;
				// ST-STDP
				float st_modulation, st_max;
				int st_type;
				float ampa_gain;
				float nmda_gain;
				float nmdaVi_gain;
				float gabab_gain;
				
				fin >> pre_percent >> min_radius_mm >> max_radius_mm >> sigma_mm >> init_noise_pct >> init_max_sum >> max >> lrate0 >>lrateN >> lstarttimems >> lendtimems >> DAstarttimems >> DAendtimems >> ampa_gain >> nmda_gain >> nmdaVi_gain >> gabab_gain >> st_type >> st_modulation >> st_max;
				//RGM: added tau_dopamine on 02/22/2012
				//cout << pre_percent << " " ;
				
				if( fin.eof() ){
					cout << "Error reading the " << Ntypes+3 << "synaptic percentages from microcircuitry.dat.  EOF reached." << endl;
					cout << "Line: " << cell_type_name << " " << layer_name << endl;
					exit(0);
				}
				
				SYN[area][type][layer][pre_area][pre_type].pct_pre_synapses = pre_percent;
				SYN[area][type][layer][pre_area][pre_type].min_radius_mm = min_radius_mm;
				SYN[area][type][layer][pre_area][pre_type].max_radius_mm = max_radius_mm;
				SYN[area][type][layer][pre_area][pre_type].sigma_mm = sigma_mm;
				SYN[area][type][layer][pre_area][pre_type].init_noise_pct = init_noise_pct;
				SYN[area][type][layer][pre_area][pre_type].init_max_sum = init_max_sum;
				SYN[area][type][layer][pre_area][pre_type].max = max;
				SYN[area][type][layer][pre_area][pre_type].lrate0=lrate0;
				SYN[area][type][layer][pre_area][pre_type].lrateN=lrateN;
				SYN[area][type][layer][pre_area][pre_type].lstarttimems=lstarttimems;
				SYN[area][type][layer][pre_area][pre_type].lendtimems=lendtimems;
				SYN[area][type][layer][pre_area][pre_type].DAstarttimems=DAstarttimems;//RGM: add DA here
				SYN[area][type][layer][pre_area][pre_type].DAendtimems=DAendtimems;//add DA here
				SYN[area][type][layer][pre_area][pre_type].ampa_gain=ampa_gain;
				SYN[area][type][layer][pre_area][pre_type].nmda_gain=nmda_gain;
				SYN[area][type][layer][pre_area][pre_type].nmdaVi_gain=nmdaVi_gain;
				SYN[area][type][layer][pre_area][pre_type].gabab_gain=gabab_gain;
				SYN[area][type][layer][pre_area][pre_type].ststdp_type = st_type;
				SYN[area][type][layer][pre_area][pre_type].ststdp_modulation = st_modulation;
				SYN[area][type][layer][pre_area][pre_type].ststdp_max = st_max;
				
				
				
				fin >> first_string;  // prepare for next iteration
			}
			
		}
		
		
		
	}

	
	void print(ofstream &fout, 	vector<vector<connection> > inSYN[MAX_AREAS][MAX_TYPES][Layers] ){
		
		for(int area = 0; area < area_data.get_N_areas(); area ++ ){
			for(int i = 0; i < Ntypes; i++){
				for(int j = 0; j < Layers; j++){
					fout << synapses[area][i][j] << " " ;
				}
				fout << endl;
			}
			fout << endl;
		}
		
		for(int area = 0; area < area_data.get_N_areas(); area ++ ){
			fout << "Area: " << area_data.get_area_name(area) << endl;
			for(int i = 0; i < Ntypes; i++){
				for(int j = 0; j < Layers; j++){
					if( inSYN[area][i][j].size() > 0 ){
						fout << "layer#: " << j << endl;
						for(int pre_area=0;pre_area < area_data.get_N_areas(); pre_area++){
							fout << "pre_area#: " << pre_area << " " << area_data.get_area_name(pre_area) << endl;
							for( int k = 0; k < Ntypes; k++){
								if( inSYN[area][i][j][pre_area][k].pct_pre_synapses > 0 ){
									fout << '{';
									fout << Names[k] << " ";
									fout << inSYN[area][i][j][pre_area][k].pct_pre_synapses << " " ;
									fout << inSYN[area][i][j][pre_area][k].min_radius_mm << " " ; 
									fout << inSYN[area][i][j][pre_area][k].max_radius_mm << " " ; 
									fout << inSYN[area][i][j][pre_area][k].sigma_mm << " " ;
									fout << inSYN[area][i][j][pre_area][k].init_noise_pct << " " ;
									fout << inSYN[area][i][j][pre_area][k].init_max_sum << " " ;
									fout << inSYN[area][i][j][pre_area][k].max << " ";
									fout << inSYN[area][i][j][pre_area][k].lrate0 << " ";
									fout << inSYN[area][i][j][pre_area][k].lrateN << " ";
									fout << inSYN[area][i][j][pre_area][k].lstarttimems << " ";
									fout << inSYN[area][i][j][pre_area][k].lendtimems << " ";
									fout << inSYN[area][i][j][pre_area][k].DAstarttimems << " "; //RGM: added on 02/22/2012
									fout << inSYN[area][i][j][pre_area][k].DAendtimems << " "; 
									
									fout << "} " ;
								}
							}
							fout << endl;
						}
					}
				}
			}
			fout << endl;
		}
		
		
		
	}

	void print(ofstream &fout){
		
		print(fout, SYN);
		
	}
	
	
private:
	
	vector<string> layerNames;
	vector<string> Names;
	
	areas area_data;
	
	int Ntypes;
	
	int	synapses[MAX_AREAS][MAX_TYPES][Layers];
	vector<vector<connection> > SYN[MAX_AREAS][MAX_TYPES][Layers];
	float cell_type_percents[MAX_AREAS][MAX_TYPES];
	
};
	
#endif
